Density of states of relativistic and nonrelativistic two-dimensional electron gases in a 

uniform magnetic and Aharonov-Bohm fields 
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We study the electronic properties of 2D electron gas (2DEG) with quadratic dispersion and 
with relativistic dispersion as in graphene in the inhomogeneous magnetic field consisting of the 
Aharonov-Bohm flux and a constant background field. The total and local density of states (LDOS) 
are obtained on the base of the analytic solutions of the Schrodinger and Dirac equations in the 
inhomogeneous magnetic field. It is shown that as it was in the situation with a pure Aharonov- 
Bohm flux, in the case of graphene there is an excess of LDOS near the vortex, while in 2DEG the 
LDOS is depleted. This results in excess of the induced by the vortex DOS in graphene and in its 
depletion in 2DEG. 
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I. INTRODUCTION 

The continuous linear energy dispersion i?(k) = 
±fi.z)ir|k| of the Dirac quasiparticle excitations when the 
homogeneous magnetic field B is applied perpendicular 
to its two-dimensional (2D) plane transforms into the 
discrete Landau levels (LLs) 
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observed in graphene. Here k is the momentum mea- 
sured from K± points, eo = \Jhv\eB j c is the relativistic 
Landau scale with vp being the Fermi velocity. The spec- 
trum ([1]) is characteristic of Dirac fermions and the break- 
through in experimental studies of graphene is caused not 
only by its fabrication^, but also by the demonstration 
of its unique electronic properties^i^ that follow from the 
unusual spectrum ([T]). 

The hallmark of this spectrum is the zero energy field 
independent lowest LL whose existence does not in fact 
depend on the homogeneity of the field 4 In general the 
inhomogeneous magnetic perturbation can be presented 
as a sum of a constant (averaged over the system) field 
and a field localized in some regions of the 2D system. A 
limiting case of the perturbation can be presented by the 
Aharonov-Bohm field which is created by an infinitely 
long and infinitesimally thin solenoid. 

The purpose of the present paper is to study the elec- 
tronic excitations in graphene in a field consisting of the 
Aharonov-Bohm fiux and a constant background mag- 
netic field. As in the first publication^ where we studied 
the Aharonov-Bohm flux only, our main goal is the in- 
vestigation of the local density of states (LDOS). We find 
that the demonstrated in Ref. rather peculiar behav- 
ior of LDOS in Dirac theory with Aharonov-Bohm field 
persists in the presence of the constant background field. 
We expect that this behavior can be observed in scanning 
tunneling spectroscopy measurements for graphene pen- 
etrated by vortices from a type-II superconductor on top 



of it. We also compare the obtained expressions with the 
corresponding results for two-dimensional electron gas 
(2DEG) with a quadratic dispersion, where the singular 
behavior of the LDOS is absent. 

In practice, such a magnetic field configuration may 
be obtained when a type-II superconductor is placed on 
top of graphene. In the previous publication, we consid- 
ered an idealized picture when the vortex is single and 
there is no impact from other Abrikosov vortices. Now 
the constant background field is supposed to mimic the 
impact of the other vortices penetrating graphene. It is 
worth to stress that devices like this, with a supercon- 
ducting film grown on top of a semiconducting hetero- 
junction (such as GaAs/AlGaAs) hosting a 2DEG, have 
in fact been fabricated twenty years ago,^'^ so it should be 
possible to fabricate the graphene based devices. While 
normally the 2DEG is buried deep in a semiconducting 
heterostructure which makes the LDOS measurements 
problematic)^ the graphene surface is open to the LDOS 
measurements. While initially the STS measurements 
were done on graphene flakes on graphite^, recently these 
measurements were carried out on exfoliated graphene 
samples deposited on a chlorinated Si02 thermal oxide 
tuning the density through the Si backgatej^^ So far all 
these measurements were done in a homogeneous mag- 
netic field and showed a single sequence of pronounced 
LL peaks corresponding to massless Dirac fermions ex- 
pected of pristine graphene. 

In a wider context, the inhomogeneous vortex-like 
field configurations arise due the topological defects in 
graphene that result in the pseudomagnetic field vortices, 
see, e.g., Refs. [TTl - [l3l Interestingly, even the nonsingular 
pseudomagnetic field configuration created by a curved 
bump on flat graphene^* results in the oscillations of the 
LDOS similar to the long-distance behavior of LDOS in- 
duced by the Abrikosov's vortex.'^ As proven experimen- 
tally, a strong localized pseudomagnetic field can be in- 
duced in graphene by a strain and the corresponding LLs 
are observed in the STS measurementsj^^ 
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Thus we hope that the combination of the vortex + 
constant background field considered in the present paper 
should be useful not only for the studies that involve a 
real magnetic field, but also for the problems that involve 
the superposition of magnetic and pseudomagnetic fields. 

The paper is organized as follows. In Sec.|ll]we intro- 
duce the model Hamiltonians and discuss the configura- 
tion of the magnetic field and the regularization of the 
Aharonov-Bohm potential used in this work. Sec. IIIII is 
devoted to the nonrelativistic case, and the relativistic 
case is discussed in detail in Sec. |IVl The structure of 
both sections is the same: we consider the solutions of the 
corresponding Schrodinger or Dirac equation, which al- 
low to write down a general representation for the LDOS 
in Sees. IIII Al and IIV Al Then, a more simple analy- 
sis of the DOS is made in Sees. HITB] and HYbI w hile 
the behavior of the LDOS is studied in Sees. IIII CI and 
IIV CI In Sec. |Vl our final results are summarized. The 
method of the calculation of the LDOS is explained in 
Appendix |X1 where as an example we firstly calculate the 
LDOS in a constant magnetic field for the nonrelativis- 
tic case. Since the problem with Aharonov-Bohm vortex 
has to be treated in the symmetric gauge, the calcula- 
tion of the LDOS in Appendix El involves the sum over 
the azimuthal quantum number, which is calculated in 
Appendix [Bl The full DOS is calculated in Appendix [Cl 
The LDOS, both in nonrelativistic and relativistic cases, 
is expressed in terms of the function calculated in Ap- 
pendix |Dj The Dirac equation in the magnetic field con- 
sisting of the Aharonov-Bohm flux and a constant back- 
ground fleld is solved in Appendix |E1 

II. MODELS AND MAIN NOTATIONS 

As in Ref. 5, we consider both nonrelativistic 
and relativistic Hamiltonians. The 2D nonrelativistic 
(Schrodinger) Hamiltonian has the standard form 

Hs = ~^^iDUDl)^ (2) 

where Dj — Vj -I- ie/hcAj, j — 1,2, with the vector 
potential A, Planck's constant ft, and the velocity of light 
c describes a spinless particle with a mass M and charge 
-e < 0. 

The Dirac quasiparticle in graphene is described by the 
Hamiltonian 

Hd = -ihvFl3{-fiDi + J2D2) + A/3, (3) 

where the matrices /3 and are deflned in terms of the 
Pauli matrices as 

^ = ^3, = (ai,C^T2). (4) 

Here C = ±1 labels two unitary inequivalent represen- 
tations of 2 X 2 gamma matrices in 2 + 1 dimension, so 
that one considers a pair of Dirac equations correspond- 
ing to two inequivalent Kj- points of graphene's Brillouin 



zone. The spin degree of freedom is not included neither 
in Eq. © nor in Eq. ([3]). In Eq. (O, is the Fermi ve- 
locity and A is the Dirac mass (or gap). An overview of 
its physical origin is given in- (see also a review^^). Here 
we only point out that the presence of a finite A allows 
one to distinguish unambiguously positive and negative 
energy solutions. 

There are numerous studies of the Dirac fermions in 
the field of a singular Aharonov-Bohm vortex (see, e.g., 
Refs. flTl - figl ) and, in particular, of this vortex and a uni- 
form magnetic field^Si^i devoted to the mathematical as- 
pects of the problem such as self-adjoint extension of the 
Dirac operator. As in the previous article to avoid the 
mathematical difficulties related to a singular nature of 
the Aharonov-Bohm potential at the origin, we consider 
a regularized potentiaP^'^"^ that depends on the dimen- 
sional parameter R: 

A(r) = ^^(r)e^, A^{r) + ^9{r - R), (5) 

where r = (r, 93, z), $0?/ is the flux of the vortex expressed 
via magnetic flux quantum of the electron $0 = hc/e 
with 77 S [0, 1[. The value rj = 1/2 corresponds to the 
Abrikosov's vortex flux. The corresponding magnetic 
fleld 

B(r) = Vx A= (^B+|^(5(r-i?)^e,. (6) 

The radius R of the flux tube determines the region r > R 
where the regularized potential coincides with the poten- 
tial of the problem with Aharonov-Bohm potential, while 
for r < R it describes a particle moving in a constant 
magnetic field. The solution of the problem is found by 
matching the solutions obtained in these regions. The 
limit i? — > can be taken at the end and allows to 
avoid the formal complications. As was shown in Ref. [23l. 
the final answer does not depend on the specific form of 
the regularizing potential provided that the profile of the 
magnetic field is nonsingular at the origin. 

We also mention recent works^iiS^ where the induced 
by the Aharonov-Bohm field charge density and current 
were studied for the massless Dirac fermions. In the first 
paper an infinitesimally thin solenoid was considered. 
The regularization by a magnetic flux tube of a small ra- 
dius R as in the present work is considered in the second 
paper It is shown that in the limit i? — > the induced 
current is a periodic function of the magnetic flux irre- 
spectively of the magnetic field distribution inside the 
flux tube and whether the region inside the fiux tube 
is forbidden or not for penetration by electrons. Also 
the value of the self-adjoint extension parameter is fixed 
by the regularization. The properties of the quasibound 
states in 2DEG with parabolic dispersion as well as Dirac 
electrons with linear dispersion in the presence of a cir- 
cular step magnetic field profile were recently studied in 
Ref. m. 



3 



III. NONRELATIVISTIC CASE 

In this section, we consider the solutions of the 
Schrodinger equation 



(7) 



in polar coordinates r = (r, ip) and using them we obtain 
the full and local DOS. These results are important not 
only for comparison with the relativistic case, but also 
because the relativistic result is constructed using the 
nonrelativistic one. 



A. Solution of the Schrodinger equation and 
general representation for the local density of states 
and its limiting rj = case 

Technically, to obtain the solutions of Eq. ([7]) in the 
regularized potential ([5]), one should solve this equation 
in two regions: r < R and r > R. Since in the first 
domain, r < R, the potential is nonsingular, only a reg- 
ular in the limit r — > solution of the radial differential 
equation is admissible. In the second domain, r > R, the 
solution contains both regular and singular in the limit 
r — >■ terms. The values of the relative weights of them 
can be found by matching radial components and their 
derivatives at r = i?. Finally, it turns out that in the 
limit i? — > only the regular solution survives and the 
wave function takes the form 



imip \m+ri\/ 2 ~-y/2 j \m+ri\ 



iy), (8) 



which also follows from the Schrodinger equation with 
a singular vortex. Here, the dimensionless variable 
y = r^/{2P) is expressed via the magnetic length I = 
{hc/eBy/"^ , L"{y) is the generalized Laguerre polyno- 
mial and the normalization constant An^m is given by 



A„ 



27r/2r(n-)- |m -I- r;| 4- 1) ' 



(9) 



The corresponding to the wave function ([S|) eigenenergy 
is equal to 



hhJr , 

= — (2n-f 1 



(10) 



where the cyclotron frequency lOc = eB/Mc, the radial 
quantum number n = 0, 1, . . ., and the azimuthal quan- 
tum number m = — oo, . . . , —1, 0,1,..., oo. In what fol- 
lows, it is convenient to express all energies of the non- 
relativistic problem in terms of the energy Eq = hiOc/2. 

Having the wave function, one can calculate the LDOS 
using the representation 



N{r,E,B) 



E 



(11) 



In contrast to the previous article;^ the presence of a 
constant magnetic field makes all energy spectra discrete, 
which demands some regularization of the S function in 
Eq. (|lip. For this purpose, we introduce a widening of 
the LLs to a Lorentzian shape: 



S{E - En^, 



ilm— 

TT En,, 



E-iT' 



(12) 



where F is the LL width. Such a simple broadening of 
LLs with a constant F was found to be a rather good 
approximation valid in not very strong magnetic fields.— 

To illustrate the method of calculation in Appendix [XI 
we derive the LDOS for the simplest case (77 = 0) of the 
constant magnetic field without vortex 



TT \ 2 nUJr 



(13) 



En = fuVc I + 2 



Here, = M/{2TTh^) is a free DOS of 2DEG per spin 
and unit area and we omitted the r dependence of the 
LDOS, because it is absent in the homogeneous field. 
One can readily obtain Eq. (jl3p in a much simplier way^ 
starting from the usual Landau spectrum 



(14) 



which follows from the spectrum (jlOp for 77 = 0, when one 
relabels n + {\m\ +m)/2 n. Here, the relabeled n cor- 
responds to the LL index rather than the radial quantum 
number. Nevertheless, in Appendix|X]we proceeded from 
Eq. ([TU]) to illustrate how to deal with a spectrum that 
is also dependent on the azimuthal quantum number m. 
As seen in Fig. [T] (a) on the dashed (red) curve, Eq. 
describes the usual quantum magnetic oscillations of the 
DOS. One can extract them analytically using the re- 
fiection formula (lAlip . Integrating the DOS over the 



- exp 



(^) 



energy with the thermal factor In 

can obtain the thermodynamic potential, whose deriva- 
tive with respect to the magnetic field gives magnetiza- 
tion. The corresponding oscillations of the magnetization 
are known as the de Haas- van Alphen effect,^- 

In a similar fashion we obtain in Appendix |^ the ex- 
pression for the LDOS perturbation, AN^{r,E,B) ^ 
N^{r, E, B) - N^{v, E) induced by the vortex 



4 



a 



/V/(r,E,B)/A/o^ 




FIG. 1: (Color online) The normalized full LDOS 
N^{r, E, B)/Nq as a function of E in the units of hojc- (a) 
ri = (no vortex and LDOS is r independent) and rj = 1/2 
for r = I. (b) Both lines are for rj — 1/2, r — 0.51, and r — 51. 
In all cases, the width is F = O.OSfioic- 



AN^{r,E,B) = - 



M sin mj 



Im 



p-ycoth{S+P) 
I _ g-2(5+/3) 

I 



du!e~^ sinh(5+/3) 



,-»;(5+/3+tj) 



1 + g-(5+^)+cj) 



il5] 



Here, N^{r,E,B) is the LDOS in the presence of the 
constant field and vortex and iVf (r, E, B) is the LDOS in 
the constant magnetic field without vortex (the argument 
r is present to distinguish the LDOS from the DOS). This 
expression has to be calculated for z > Q with the analytic 
continuation z — )■ —{E + iT)/Eo done at the end of the 
calculation. The representation (|15p for the LDOS is our 
starting point for the analysis of the LDOS and DOS. 
Next, in Sec. IIII Bl we begin with a simpler case of the 
DOS and return to the LDOS in Sec. UlTC] 



2D volume (area) of the system factor I/2D, this is not so 
in the presence of the vortex when the LDOS is position 
dependent. Then the full DOS per spin projection is ob- 
tained from the LDOS pT|) by integrating over the space 
coordinates 



N^iE,B) 



2ir 



dip / rdrN,j(r,E,B). 



(16) 



B. The density of states 

The details of the derivation of the DOS difference. 
While in the constant magnetic field the LDOS is po- AN^{E, B) = N^{E, B) - N^{E, B) with N^{E) being 
sition independent and is related to the full DOS by the the full DOS in the presence of the constant field without 
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vortex, are given in Appendix [Cj We obtain 



where 



X -0 



1 



1 E + iV 

2 huJr 



-Im 



E + iT 

hide 

E + iT 



V 



(17) 



Since the digamnia function ip{z) has simple poles for 
z ~ 0,-1,-2,..., it is easy to see in the clean limit 
r — 0, that the DOS difference (|17l) reduces to a set of 
5 peaks corresponding to the LLs: 



AN^{E,B) 



oo 



1 — ri)5 ^E — huc ( n 
1)5 \E-niu. 



n=0 



(18) 



The physical meaning of (ITSl) is that^^ on each LL £"„ — 
hijjc{n + 1 /2) , n + l — rj states disappear and n + 1 appear 
at the energy En = hajc{n + 1/2 + 77). 

The limit of zero field, i? — ?> 0, can be obtained from 
Eq. (jl7p using the asymptotic expansion 



ihiz) = In Z ;7 

' 2z 12z2 



(19) 



Then in the limit F — >■ 0, we reproduce the Aharonov- 
Bohm depletion of the DOS^i^Si^a at the bottom of the 
spectrum 

AN^{E, B = Q) = N^{E, B ^ 0) - V2dN^ = 

1 (20) 

caused by an isolated vortex. Integrating Eqs. (PH)) and 
(lisp (with an appropriate regularization) one can check 
that the total deficit of the states induced by the vortex 



dEAN^{E,B) 



-77(1 - r,) (21) 



does not depend on the strength B of the nonsingular 
background field. 



C. The local density of states 

The regularization parameter 6 in Eq. (|15p is important 
for the calculation of the DOS made in Appendix iBl the 
integrand of Eq. (jlSp remains regular even in the limit 
(5—^0. Therefore we can take this limit and rewrite 
Eq. (fT5)) as follows: 



AN^{r,E,B) 
X Im 



M sin vrr/ 
(7rft)2 27r 



E + tT 

I {y,z^ — E — 



(22) 



poo ^— ycoth/3 



(23) 



and the variable y describes the spatial dependence. Al- 
though the integrals in Eq. can be evaluated numer- 
ically, this computation becomes troublesome when the 
analytic continuation from z > to the complex values 
z — > —{E + iT)/ Eq is done before the numerical integra- 
tion. Thus our purpose is to derive such a representation 
for I{y, z, ri) that it can be easily computed after the ana- 
lytic continuation is done. The function I{y, z, 77) is found 
in the Appendix |D] and is given by 



r ( ^(z + 2i^-i 



^ i^(l-z-r,)/2,(l-r,)/2(2/) 



+r ('^"l r ( " + F(2_,_,)/2,,/2(j;), 



(24) 



where the function Fx^i^^ijj) is given by Eq. (IDlSp . 

The results of the numerical computation of the LDOS 
on the basis of Eqs. and (IMl) are shown in Figs. [1] and 
[2] We emphasize that in Fig. [1] we plot the full LDOS 
N^{r, E, B) as a function of energy E for fixed values 
of r, and in Fig. [21 the same quantity is presented as a 
function of the distance r from the vortex center for fixed 
values of E. Since Eq. (1^^ describes the perturbation 
of the LDOS AN^{r, E, B) by the vortex, to obtain the 
value of the fuh LDOS N^{r,E,B), we add to AN^ its 
7^ = value, which is given by Eq. (fTS]). We note that in 

A/i/2^(r,E,B)/A/o^ 




FIG. 2: (Color online) The normalized full LDOS 
Ni/2{r, E, B)/No as a function of the distance r measured 
in the units of the magnetic length I for four values of 
E/hhJc — 0.5, 1.5 (usual LLs) and E/fkUc = 1,2 (vortex-like 
levels). The width is L = 0.05/kJc. 

contrast to Ref. 5, when plotting these figures, we did not 
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take into account the presence of the finite carried density 
in 2DEG by shifting the energy origin. This makes more 
straightforward a comparison with the Dirac case, where 
low carried densities are indeed accessible experimentally. 

Although the model we consider is suitable for all val- 
ues of the distance from the center of the vortex r, there 
are obvious physical limitations on the possible value of 
r if the vortex penetrating graphene is coming from a 
type-II superconductor. First of all, r cannot be smaller 
than the vortex core, which is at least on the order of 
magnitude larger than the distance scale ro of the order 
of the lattice constant. We remind that in the previous 
paper ^ the distance r was measured in the units of rg, 
because for B = there is no such natural scale as a 
magnetic length. Secondly, we replace the magnetic field 
created by the other vortices by a constant background 
magnetic field. This approximation may be appropriate 
if one considers a vicinity of the selected vortex, which 
implies that r has to be less than the intervortex dis- 
tance ly. This distance is proportional to the magnetic 
length j-'^- Iv = Cy/irl w 1.77Z, where c « 1 is the geomet- 
ric factor dependent on the Abrikosov's lattice structure. 
Thus although one can investigate the regime r ^ I the- 
oretically, in practice it is not accessible. 

In Fig. [T] (a) we compare the already discussed after 
Eq. ([T5)) case of the constant magnetic field with the case 
when the Abrikosov vortex is also present (77 = 1/2) for 
r — I. We observe that while for 77 = [the dashed 
(red) curve is, obviously, r- independent] only the peaks 
at half-integers E/huc are present, for t] = 1/2 the weight 
of these peaks is reduced and a set of the new peaks at 
the integers E/hujc on the solid (blue) curve is devel- 
oped. This behavior can be foreseen from the expression 
for t he ful l DOS difference (HH) [or Eq. (H?])] discussed in 
Sec. IIII CI The case with the Abrikosov vortex is further 
explored in Fig. [1] b, where we plot the energy depen- 
dence of the LDOS for r = 0.51 [the solid (blue) curve] 
and r = 51 [the dashed (red) curve]. Comparing the 
results for r/l = 0.5, 1., and 5.0 we find that as the dis- 
tance r decreases, the integer E/fuvc peaks are getting 
stronger, while for r = 5.01 they practically disappear. 
This behavior allows to attribute the corresponding en- 
ergy levels to the vortex. On the other hand, the half- 
integer E/hjJc peaks corresponding to the usual LLs 
formed in a constant magnetic field are getting weaker 
as the distance r decreases. We stress that even for an 
arbitrary vortex flux 77, the latter levels will not change 
the positions, while the levels related to the vortex will 
shift their energies. 

Analyzing Eq. ([M)) . which was used to plot Fig. [TJ we 
observe that the positions of all peaks are controlled by 
the gamma functions r(z) which contain simple poles 
for z = 0,-1,-2,.... However, the intensity of the 
peaks depends on the rather complicated modulating 
function F\^fj^{y). For example, we verified that despite 
that the gamma function T[{z — l)/2] in the second term 
of Eq. ([M)) contains the pole at the negative energy 
E = —hujc/1, the final LDOS does not contain this pole. 



To gain more insight on the behavior of the LDOS we 
have investigated its behavior in the limits r — >■ and 
r ^ CO. Taking into account the y — >■ limit of Im/ given 
by Eq. ((D23| . we obtain that the value AiV^(r = 0, E, B) 
is equal to the negative LDOS p3)) in the constant mag- 
netic field. This implies that the full LDOS in the center 
of the vortex is completely depleted. 



N^{r^O,E,B) 



0. 



(25) 



Formally, this property refiects a simple fact that all solu- 
tions (IS]) of the Schrodinger equation are vanishing at the 
origin, ipn.mir = 0,(^5) — 0. This vortex-induced deple- 
tion of the LDOS in the nonrelativistic 2DEG was already 
seen in Ref. and now we conclude that it should also 
occur in the presence of the background magnetic field. 
This is exactly what we observe in Fig. [21 where all four 
curves begin from zero. Two of these curves, viz. the 
solid (blue) and the dash-dotted (black) are for the usual 
LLs with E/huc ~ 0.5,1.5, and the other two [dashed 
(red) and dotted (violet)] are for the vortex levels with 
E/hujc = 1,2. For small r < / all curves increase lin- 
early as expected from the analytic results described in 
Appendix [Cl if we take there 77 = 1/2. Since for the large 
y the function Fx^ decays exponentially [see Eq. (jD26|) ]. 
the LDOS difference AN^{r, E, B) - e-'''/^'" for r -> 00. 
Accordingly, the large r behavior of the full LDOS de- 
pends on the contribution of the position independent 
LDOS (US)). Thus the large r limit of ah curves in Fig. [H 
is determined by the corresponding value of the LDOS 
in the dashed (red) curve in Fig. [T] (a). 



IV. RELATIVISTIC CASE 



In Sec. IIV Al we consider the solutions of the Dirac 
equation 



where the wave function is now a spinor 



*(r,C) = 



V'i(r,C) 
V'2(r,C) 



(26) 



(27) 



and the index C labels two inequivalent K-j- points. No- 
tice that in Appendix |E1 the definition (jEll) for '4>2 ex- 
plicitly includes the factor i. Using these solutions in 
Sec. IIVBI we obtain the full and the local DOSs that is 
considered in Sec. IIV Cl 



A. Solutions of the Dirac equation and general 
representation for the local density of states and its 
limiting 77 = case 

The Dirac equation (j26|) with the regularized potential 
^ is solved in Appendix |E1 A general strategy is the 
same as described in Sec. IIII Al but the main difference is 
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in the matching conditions. While the radial components 
of the spinor ^>{r) have to be continuous: 



MR + OX) ^ MR- 0,0, 
MR + o,0 = MR-o,0, 



(28) 



their derivatives in contrast to the nonrelativistic case 
have a discontinuity: 



MR + 0, - i^'iiR -0,0 = ^MR, 0, 
MR + 0,0- ^2{R -0,0 = -^MR, 0- 



(29) 



The discontinuity of the conditions (1^^ follows from 
Eq. (jE2l) with the discontinuous potential (l5|). Another 
way to apprehend this discontinuity is to obtain a sin- 
gular SLT r = R pseudo-Zeeman term squaring the Dirac 
equation (see Ref. Isj for an overview) . 

After the limit ii — )■ is taken, we obtain the following 
solutions: 



fr 1) = 



1 



2^ -\/ T^Sn.r. 

for m > 0, 



^£n, m±Ae'("^ -^^^J:^+,,My) 
±z^£:„,„TAe™'^j;;,+^(y) 



(30) 



*i'b^(r,l) = 



1 



for m = 0, and 
^f^^^) fr 1) = 

1 



Vg.„, o±Ae-^ '^Ji%(t/) 
T»V^„,oTAJ!!+i(y) 



(31) 



2ly^TT£„ 



Vg„. „±Ae^("' -i)'^J|"^^^_,|(y) 
T*V^„,™TAe™'^J|"J;i^|(y) 



(32) 



for m < 0. Here, the upper and lower signs ± correspond 
to the positive and negative energy solutions, 5^^^ — 
±£n,m with the absolute value of the energy 



An,. 



A2 + egA,, 



2n + \m + rj ~ 1\ + m + rj + 1 



(33) 



for n > 0, the function J"(j/) is defined by Eq. (|E14[) 
with y as in the nonrelativistic case, and the relativistic 
Landau scale eo is defined after Eq. ([1]). The zero- mode 
solution with £ = — A is a holelike 



<™(r,l) = 



1 



2nl 







Let us now compare the solutions of the Schrodinger and 
Dirac equations with the zero azimuthal number, m — 0. 
One can see from Eq. (jS]) that ^nfiif, 1) ~ r^, because^ 



KiO) 



n + a 
n 



r{a + n + 1) 
r(a + l)n! 



(35) 



On the other hand, from Eqs. (|3T|) and p4l) for m = we 
observe that while the upper components are regular at 
r = 0, the lower components diverge as ^2n,o{''', 1) ^ r^^K 
Comparing these results with the behavior of the wave 
function in the Aharonov-Bohm field'^ we observe that 
the presence of the background magnetic field does not 
change the asymptotics of the m — solutions for r — > 0. 
Also as expected,^ the zero-mode solution ([M)) for ( ~ 1 
and chosen direction of the field is holelike, £ = —A. 
The solutions for the case ( = —1 are the following: 



^f^^^' fr -1) = 

1 



2ly/7T£n. 

for m > 0, 



Ty/£r 

Vfn,mTAe*(™-i)'^J,".+,-i(2/) 



-'n,m i A e"""^ J^.|_^ (y) 



(36) 



<'b^(r,-l) 



1 



2ly/n£nfi 
for m — 0, and 



±^£n.o± AJ^+\y) 
i^£n,oTAe-^'^JMiy) 



(37) 



^ n,m V* , 
1 



-1) 



2M-K£n,'. 



±Vg„,^± Ae^-'^J|"J;^^|(y) 
zVf„,™TAe^(™-i)'^J|';„^^_^|(y) 



(38) 



for m < 0. Again the signs ± correspond to the solutions 
= ±£„^„i with 71 > and the energy given by 
Eq. ([55)) . Now the zero-mode solution 



«'o,m(r,-l) 



1 



2ttI 



„imip jQ 





(y) 



m < 0, (39) 



is electron-like, £ = A. Also the lower components of the 
m — solutions ([57)1 and ([5^ are regular at r = 0, and 
the upper components diverge as ipin,o{f, —^) ^ r^^- 

Since the solutions of the Dirac equation are character- 
ized not only by the quantum numbers, but also by the 
sublattice label A and B, energy ±, and the valley index 
C = ±1, instead of directly writing an analog of Eq. (|lip . 
it is more convenient to construct the Green's function 
expressing the LDOS via the combinations of its matrix 
elements. The eigenfunction expansion for the retarded 
Green's function reads 



oo oo 



m<0. (34) G^ir,r',E + iO;0=Y. E 



*i+l(r,C)*ltm(r',C) 



E - £n,m. + iO 

^ U(i-,C)^d/(rO \ 

E + £n,m +iO j ' 



n— m— — oo 



(40) 



The LDOS for A and B sublattices is expressed in terms 
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of the Green's function (l40l) as follows 



--Im[G„ii(r,r,i; + ir;C = 1) 

TT 

+ G^ii(r,r,i; + zr;C- -1)] 
--Im[G„22(r,r,i; + ir;C = 1) 



(41) 



-1)] 



where similarly to the nonrelativistic case, the LL width 
r is introduced. Substituting the solutions of the Dirac 
equation in the Green's function (UHl), and using the def- 
inition (|4T|) . we obtain 



1. 



X G[y,z- 



E + ir±A 



(E + ir)2 - A2 



(42) 



where the upper (lower) sign corresponds to A (B) sub- 
lattice, the relativistic Landau scale eo is defined be- 
low Eq. ([T|), and the normalization constant = 
eg/ {2Trh'^Vp) corresponds to the value of the free (ry = 
B = T = 0) DOS for the Dirac quasiparticles per spin 
and one sublattice (or valley) taken at the energy E — 
when the energy gap A = 0. Writing Eq. we intro- 
duced the function 



3 

E 



(43) 



which consists of the three terms 



g2{y,z,'i]) = 



oo oo 

E E 

11—0 m— — oo 

oo oo 



\m+ri\ ^ 
n— u m— — oo 



(44) 



E 

n=0 



z + 2{n + rj) z + 2n 



with \n,m defined in Eq. ([33]) and 
K.,n = 2n + \m + 'nl - 



(45) 



Note that the function 173 contains the singular terms that 
originate from m = solutions of the Dirac equations. 

The gi,2 contributions are calculated in the same way 
as was derived Eq. ()15p in Appendix [Aj viz., exponen- 
tiating the denominators [see Eq. (|A2|) ] and introducing 
the regularizing factor S > 0, and then using the sum 
^41). we obtain 



gi{y,z,ii) = - / 

Jo 



d/3e 



-f}z_ 



I _ g-2(5+/^) 

V p-('5+/3)(™+'j);-, 



-acoth(5+/3) 



sinh(^ + 



(46) 



where we also shifted the dummy index m — )■ m + 1. 
As we saw in the nonrelativistic case, the presence of d 
is necessary for the calculation of the DOS, although it 
can be omitted in the expressions for the LDOS. The 
remaining sum over m can be found using Eq. (|B7[) from 
Appendix |B] 



Agi(2/, z, r/) = gi{y, z, 77) - gi{y, z, 0) 



Sm TTT] 



-f}z_ 



-2{S+0) 



1 



-2(<5+/3) 



-ycoth{5+P) 



(47) 



duje^^ cosh uj/ sinh((5+/3) 



1 + g-{u:+S+l3) 



where we introduced the function Agi, which describes 
the perturbation by the vortex. Similarly, for Ag2 we 
have 



^92{y, z, 77) EE 52(2/, z, 77) - 32(2/, z, 0) 

^sin^ 'e~ 

■K Jo i_e-2(<5+/3) 



ycoth(5+l3) 



(48) 



duje~^ cosh sinh((5+;3) 



The case of gs is even simpler, because there is no sum- 
mation over 771. Using the sum (jA4p we obtain an analog 
of Eq. pS)) . It contains the difi^erence of two modified 
Bessel functions, which can be expressed via the Mac- 
Donald function™ 



2 sin TTi' 

Finally we arrive at the result 
2 sin 7r77 



[/_,(x)-/,(a;)]. 



(49) 



93{y,z,T]) = 

/•oo 

X / d/3e 
Jo 



-I3z_ 



-2{5+l3) 



g-ycoth(5+/3)^ 



" Vsinh((5-f-/3) 
(50) 



Notice that since g3{y,z,ri = 0) = 0, there is no need 
to introduce a function Ag^. Having the functions A5'i_2 
and 33 we can directly calculate the LDOS perturbation 
by the vortex, AN^^'^'^\r, E, B) = N^^^'^\r, E, B) - 



N, 



DiA,B] 







{r,E,B). 



The subsequent consideration is made in parallel to 
the nonrelativistic case. We consider first the LDOS in 
the constant magnetic field {rj = 0) when due to the 
translational invariance it coincides with the DOS per 
unit area. The LDOS can be derived in a similar to 
Eq. (IT^ way, but a special care has to be taken because in 
contrast to the nonrelativistic case, the cutoff parameter 
6 enters the final result. 



n^^^-^\e,b) 



X Im 



E + iT±A 



£0 



ln(2<5) + 7 + (I) + ^ 



(51) 
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where we kept only the divergent in the hmit 5^0 
terms It is more convenient to rewrite Eq. (jSip in the 
form derived in Ref. [IM where instead of the cutoff 5 the 
bandwidth W cutoff is used 



N, 



B{A,B) 







2 



r , 

— in TT^ 







A^-{E + iV) 



Im 



4 



E + ir±A 



(52) 



The advantage of the representation (|52l) is that its B — 
hmit takes the usual form.-* The quantum magnetic 
oscillations of the LDOS, N^{E,B) ^ N^'''^\e,B) = 
iV°^^^ {E, B) for a = are shown in Fig. H (a) on 
a dashed (red) curve. Only the positive-energy region 
is shown, where the positions of the peaks, £n/^o = 
^/2n, are in accord with the Dirac spectrum (H]). The 
nonequidistant LLs along with the peak at = related 
to the energy independent lowest LL are characteristic of 
the Dirac fermions. The reflection formula (|A11I) allows 
one to extract these oscillations analytically 32, 

The representation (|^^ for the LDOS, where the func- 
tion ga]) consists of the three terms gT]), (gH]), and ((50)) 
which describe the LDOS perturbation, is our starting 
point for the analysis of the LDOS and DOS in the rel- 
ativistic case. In the next Sec. IIVBI we begin with the 
DOS and in Sec. |IVC] return to the LDOS. 



The density of states 



A/,°(r,E,B)/A/o° 

7' 




A/i/2°(r,E,B)/A/o° 




FIG. 3: (Color online) The normalized full LDOS 
N^{r, E, i?)/A^°(eo) as a function of energy E in the units of 
the relativistic Landau scale eo- The LDOS is an even func- 
tion of E, so only the positive-energy region is shown, (a) 
?7 = (no vortex and LDOS is r independent) and ri — 1/2 
for r — I. (b) Both lines are for r; = 1/2, r = 0.51, and r — 41. 
In all cases, the width is F = 0.05eo, W/eo — 3.35, and A = 0. 



The full DOS per spin projection is given by the spa- 
tial integral ([T51) . Accordingly, the full DOS perturba- 



tion by the vortex AN^^'^'^\e,B) 



N, 



D{A,B} 



D(A,B] 



{E,B) 







{E, B) for A and B sublattices takes the form 



E + iV±A 



(53) 



dy{Agi{y,z,r]) + Ag2{y, z,t]) + gsiy, z,r])) 



The subsequent calculation on the basis of Eq. ([531 
is similar to the nonrelativistic case considered in Ap 
pendix[Cl and gives [compare Eqs. (|C6|) and (fTTl) ] 



+(z + 2r,) ("^(1) 



— Im 



E + iT±A 

2^^1 
z + 2rj 



2,(1 + 1 



(54) 



where z -[{E+ir)'^-A'^]/el. In the clean limit F -> 
the DOS difference reduces to 

AN^<-'^-^\e,B) ^7^diE±A) 



■ 2{E ± A)sgn£; 



J2 nSiE^ - A2 - 2(71 + Tj)el) 



n=l 



(55) 



where except of the first, proportional to ry, zero-mode 
term each S function corresponds to both positive and 
negative energy peaks. A comparison of this result with 
Eq. (|18l) for the nonrelativistic problem sheds the light 
on the difference between these cases. We observed from 
Eq. (|T8)l that all peaks associated with the usual LLs are 
depleted, while the peaks related to the vortex are de- 
veloped. At first sight, Eq. ([55]) follows the same pat- 



tern, viz. the LL peaks with = zL^/A^ 

with n — 1,2,... are depleted and the vortex-like lev- 



2eln 
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els S^"^^ = ±v/A2 + 2e2(n + r;) with n = 1, 2, 3, . . . are 
developed. However, the first term r]S{E ± A) related to 
the zero-mode solutions of the Dirac equation is present 
for any magnetic field configuration and the addition of 
the vortex only adds rj to the weight of the corresponding 
peak. This property is an illustration of the topological 
origin of the lowest LL.'^^ 

The _B — limit can again be obtained using the 
asymptotic expansion ()19|) . which for the expression in 
the square brackets of Eq. ([M)) gives 



277 1 



+ (z + 277) (V' (I) - V' 



z + 2r] 



2rf 

z 



O 



(56) 



Substituting Eq. (|56p in Eq. and making the analytic 
continuation z — > —[{E + iT)'^ — A^]/eg in the clean limit 
r — >■ 0, we obtain 



D{A,B) 




AN^(^'^\e,B = 0) = N^^^'^\E,B = 0) - V2dN, 
= rf5{ET^)- 

(57) 

This result is in agreement with Refs. |3^ and where 
the DOS p^{E, C,) for a separate K± point, but summed 
contributions for A and B sublattices, was considered. 
Its perturbation Ap^{E, () = p^{E, () - V2dPo{E) with 
respect to the free DOS per spin and one valley, p^{E) — 
\E\9{E^/A^ - l)/{2TTfi^vj,), is equal to 

Ap^iE, = - ^77(1 - tjME - A) + SiE + A)] ^^^^ 

+ 77^(£; + CA), 77 >0. 

Integrating Eq. (157)) one can find the total excess of 
the states induced by the vortex 



AA^, 



D _ 



/oo 
dE{AN^^'^\E,B)+AN^^'^\E,B)) = 2rf 
-OO 



(59) 



As in the nonrelativistic case (PT|) . it turns out that the 
integral ([5^ does not depend on the strength B of the 
background field. This can be checked by integrating 
the sum (j55p and using an appropriate regularization. 
Completing our discussion of the DOS we note that the 
value AN^ has to be distinguished from the induced by 
the magnetic flux fractional fermion number— which in 
terms of the DOS (l58l) can be written as follows: 



1 Cn 
-- J^^dEsgnEAp^{E,0^f. 



(60) 



function /(y, z, rf) defined by Eq. (|23|) . which was used in 
Sec, nil CI to express the nonrelativistic LDOS 



Agi(y,z,'q)^ /(7/,z + 1,77), (61) 



and 



Ag2{y,z,r/) 



sm vrr/ 
2tt 



iiy,z~- 1,77). 



(62) 



Thus the only remaining term we have to find is given 
by Eq. ([5(7|) . Changing the variable x = , we obtain 



93{y,z,v) 



2 sin TTT] 



1 — a;^ 



2xy 

1- x^^ 



(63) 



Now using the integral (2.16.10.5) from'^'^ (one can also 
change the variable to t via — [t/{l + 1)]^/^ and use 
the integral (jD9[) ). we can write 



dx 



J.2 _ y2 



exp I -b— I K„ 

y^ — X 



2cx 

y2 _ J.2 



4c \ 2 J \ 2 

X W^(i-a)/2,./2 (b+^b^-ic/yy) 

X W^(l-a)/2W2 (b~Vb^-{c/yr 



(64) 



Thus we can express 53 in terms of the Whittaker func- 
tion W\fj,{z) as follows: 

gaiy, z, v) = ^^ I^{y, z, ry) (65) 



2tt 



with the function 
1 



(66) 

Thus the final expression for the relativistic LDOS per- 
turbation by the vortex takes the form 



E + iT±A 



X AG[y,z^ _(^+^r) A_^^ 



(67) 



where the function 

AG(,,z,77) = ^ (68) 
X [Hy, z -I- 1, 77) 4- I{y, z - 1, 77) + I^{y, z, rj)] 



C. The local density of states 

The contributions A51 2 to the relativistic LDOS given 
by Eqs. (H7)) and can be written in terms of the 



is expressed via the defined above functions p4|) and ([66]) . 

To complete the analytic treatment, we consider the 
behavior of the LDOS in the most interesting case of the 
small r, when we expect that the difference between the 
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relativistic and nonrelativistic cases should be the most 
transparent. The observation (|25p that in the nonrela- 
tivistic case the full LDOS in the center of the vortex 
vanishes turns out to be useful for better understanding 
of the relativistic case. Indeed, let's consider the first 
two terms of Eq. (|43l) with gi^2 that contribute to the 
fuU LDOS (gl]). Since the numerators of 51,2 in Eq. pi)) 
vanish at y = 0, the only term that governs the behav- 
ior of the full LDOS in the r limit is the function 
I^{y, z, r]), which due to its origin from the m — solu- 
tions is expected to be divergent. 

The same result can be verified using the final ex- 
pressions §31 and §^ for A7V°''^'^^(r = 0,E,B). For 
y = 0, the first two terms of Eq. ([68|) with the function /, 
which originate from Agi^2 [see Eqs. (pT|) and (15^ ] can 
be combined together; 



Sm TTT] 

— — Im 



Im[/(?; = 0, z + 1, r/) + I(y = 0, z - 1, ??)] 



(69) 



where we used the value Im/(y = 0,z,i]) established in 
Eq. (jD23l) and then transformed the first digamma func- 
tion using Eq. (jD22[) . Thus we find that in the limit 
r — >■ the contribution of these Agi,2 terms to the LDOS 

difference AA^^^^^'^^r = 0,E,B) given by Eq. §7^1 is 
equal to the negative LDOS ((5T|) in the constant mag- 
netic field. 

Let us now analyze the behavior of the function 
I^{y,z,ri) in the r limit. Using the expansion of 
the Whittaker function (jD19p in the limit ?/ — ?> 0, we 
obtain 



_r(./2r(,) 

r(r; + z/2) ^ yy h y 

Thus the full LDOS is divergent at the origin as 



(70) 



E + ir±A 



?f2 



9f2 



with the case when the vortex is also present (77 = 1/2) 
for r — I. Since we consider the situation when A = 0, 
there is no difference between sublattices, Nq{E,B) = 



D{A) 



{E,B) = Nq'^^\e,B) and the LDOS is an even 
function of energy, so the positive energy region is plot- 
ted. We observe that compared to 77 = [the dashed 
(red) curve] for 77 = 1/2 [solid (blue) curve] a set of the 
new peaks at £„/eo = \/2n + 1 with n = 1, 2, ... is de- 
veloped and the lowest LL peak (n = 0) is enhanced. 
This behavior can be foreseen from the expression for 
the full DOS difference Eq. ([35]) [or Eq. ^] discussed 
in Sec. IIVBI The case with the Abrikosov vortex is fur- 
ther explored in Fig. [3] (b), where we plot the energy 
dependence of the LDOS for r — 0.51 [the solid (blue) 
curve] and r = Al [the dashed (red) curve]. Comparing 
the results for r/l = 0.5, 1.0, and 4.0 we find that as the 
distance r decreases, the peaks at £n/^o — \/2n -I- 1 with 
n = 1,2,3... related to the vortex are getting stronger. 
When r further decreases, the peaks related to LLs grow 
faster than the vortexlike peaks. This behavior indeed al- 
lows to attribute the corresponding energy levels to the 
vortex. On the other hand, the peaks £n/^o = \/2n with 
n — 1,2,... corresponding to the usual LLs (H]) are get- 
ting weaker as the distance r decreases. We remind that 
even for an arbitrary vortex fiux 77 the latter levels will 
not change the positions, while the levels related to the 
vortex will shift their energies. Fig. [3] also illustrates a 
special character of the lowest LL that is present even in 
an inhomogeneous magnetic field (see also recent simu- 
lations in Ref. IstI ) , and therefore is getting stronger as r 
decreases. 



(71) 




For 77 — 1/2, the divergence is 



as was in the ab- 



sence of the background fieldi^ As we discuss below, the 
presence of this field makes the divergence of the LDOS 
strongly energy dependent. 

The results of the numerical computations of the full 
LDOS on the basis of Eqs. (1571) and ([55)1 are shown in 
Figs. [31131 andjSl Since Eq. (1571) describes the perturba- 
tion of the LDOS A7V°(r, E, B) by the vortex, to obtain 
the value of the full LDOS 7V°(r, E, B), we add to AiV,^ 
its r? = given by Eq. ([5^ . 

In Fig. 131 (a) we compare the already discussed after 
Eq. (15^ case of the LDOS for a constant magnetic field 



FIG. 4: (Color online) The normalized full LDOS 
N^l'^'^\r, E, B) /N^ {eo) as a function of energy E in the 
units of the relativistic Landau scale eo for r — I. The gap is 
A = eo, the width is T = O.OSeo, and W/eo = 3.35. 

In Fig. |4l we consider the energy dependence of the 
LDOS Nf^^^'^\r,E,B) when there is a gap A = eo 
in the spectrum. The distance from the vortex center 
is r = L The gap introduces asymmetry between the 
LDOS on A and B sublattices and also makes the LDOS 
asymmetric with respect to i? = 0, so we have to plot 
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both negative- and positive-energy regions. Indeed we 
observe that the zero LL peak at £ = A is present only 

in N^I'^\r,E,B), while the peak dX E — —A shows up 

only in Nfl'^\r, E, B). The vortexlike levels also be- 
come asymmetric with respect to i? = 0. All this illus- 
trates that the STS on graphene on a substrate that can 
induce inequivalence of sublattices in graphene should 
reveal these features. 

From Eq. ([7T|) we expect that the presence of the back- 
ground magnetic field makes divergence at r — 
of the LDOS strongly energy dependent: it is empha- 
sized by the poles of the first F function, when the 
energy E is close to the energies of the usual LLs, 

and oppo- 



difference AN^{r,E, B) 



, the large-r behavior 



iy^A^ 2eln with n = 0,1,2,.. 



sitely, because the second F function is in the denomi- 
nator, when E is equal to the energies of the vortexlike 
levels, fi^^ = ±^A2 Jr2el{n + ri) with n = 1, 2, . . . , the 
divergence is suppressed. This is exactly what we ob- 



A/i/2"(r,E,B)/A/o 




FIG. 5: (Color online) The normalized full LDOS 
NF/2(J'j B)/No' {eo) as a function of distance r measured 
in the units of the magnetic length I from the vortex for 
four values of E/eo = V2, 2 (usual LLs) and E/eo = v^, V5 
(vortex-like levels). The width is F = 0.05eo, W/eo = 3.35, 
and A = 0. 

serve in Fig. [SI where we show the dependence of the 
LDOS on the distance r for fixed values of the energy 
(A = 0). Indeed, the solid (blue) and dash-dotted (black) 
curves which correspond to energies E/eo = V^, 2 of the 
usual LLs have divergent behavior at the origin. Ob- 
viously, this divergence is also present for E = and 
the corresponding curve will be above the higher energy 
curves E/eo = and E/eo = 2. On the other hand, 
the dashed (red) and dotted (violet) curves, which cor- 
respond to the energies E/eo — V^, and ^/b of the vor- 
texlike levels tend to go to a constant value at r = 0. 
Strictly speaking divergence is suppressed only when 
the function F~^(77 4- z/2) in Eq. ([7T|) is zero, but for 
small values of the level width F, the divergence seems 
to be completely suppressed for the chosen values of the 
energy. For r ^ I, the behavior of the LDOS resembles 
the nonrelativistic case. Since in this limit the LDOS 



of the full DOS is determined by the contribution of the 
position independent LDOS ([5^. Thus the large r-limit 
of all curves in Fig. [5] is determined by the corresponding 
value of the LDOS in the dashed (red) curve in Fig. [3] (a) . 



V. CONCLUSIONS 

The main motivation of this work was to address the 
question as to whether one can distinguish graphene from 
2DEG by measuring the LDOS near the Abrikosov vortex 
penetrating them. In the first publication,^ we investi- 
gated the simplest formulation of the problem with a sin- 
gle vortex. In the 2DEG, the solutions of the Schrodinger 
equation in the presence of the Aharonov-Bohm field are 
regular and the LDOS near the vortex is depleted. On 
the other hand, a specific feature of the Dirac fermions 
in the field of the Aharonov-Bohm flux, namely, such as 
the presence of the divergent as r^^ at the origin the 
m = solution of the Dirac equation, results in the r~^^ 
divergence of the LDOS in the vicinity of the vortex. 
Therefore the LDOS enhancement near the vortex can 
really distinguish graphene from 2DEG. 

This positive answer obtained in the previous paper- 
is now extended for the case of a more complicated 
magnetic field configuration consisting of the Aharonov- 
Bohm flux and a constant background fleld, as one can 
see just from a comparison of Figs. [Hand [5l It turns out 
that the character of the divergence in the Dirac case re- 
mains the same, but it is now strongly modulated by the 
energy-dependent factor. The divergence is present when 
the energy is equal to the energies of the usual LLs ([IJ , 
including the lowest zero-energy LL. 

The signiflcant difference between the relativistic and 
nonrelativistic cases can be understood by comparing the 
squared Dirac equation with the Schrodinger equation. 
While the Schrodinger equation contains only an effec- 
tive centrifugal potential, which originates from the an- 
gular part of the Hamiltonian, an equation for one of 
the components of the Dirac spinor always contains an 
attractive pseudo-Zeeman term. We call it the pseudo- 
Zeeman term because it is related to the sublattice rather 
than to the spin degree of freedom. Since for the zero az- 
imuthal number m the centrifugal part of the potential 
is the smallest, the attraction term results in the diver- 
gence of the LDOS near the vortex. Our main results, 
which allow to conclude that this picture remains valid 
in the constant background fleld can be summarized as 
follows. 

(i) We obtained analytic expression for the LDOS per- 
turbation by the Aharonov-Bohm flux in the presence of 
a constant background magnetic fleld in the nonrelativis- 
tic, see Eqs. ([^ and ([M)) . and relativistic, see Eqs. (p7)) 
and (j68p cases. The nonrelativistic answer is written in 
terms of the function (j24p which is expressed as a com- 
bination of the Whittaker functions in Eq. (jDlSp . The 
relativistic answer (1681) is expressed in terms of the same 
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function (|24|) . and a function ((66|) that describes the con- 
tribution of the TO = solutions of the Dirac equation. 

(ii) We show that in the vicinity of the vortex (r < 0.21) 
the relativistic LDOS is governed by the function 
so that in the hmit r — > the LDOS is given by Eq. ([7T|). 

(3) We obtained compact analytic expressions for the 
DOS perturbation by the Aharonov-Bohm flux. For the 
nonrelativistic case, this is Eq. (jl7p . which in the clean 
limit reduces to the known^^ result given by Eq. (ITSt . For 
the relativistic case the corresponding expressions for the 
DOS are Eqs. ^ and ([SS]). 

We hope that the obtained results will be useful both 
for the experimental STS studies of graphene and for the 
theoretical studies of the interaction effects in an inho- 
mogeneous magnetic field similar to the recent work.— 
Among possible extensions of the considered problem we 
mention the necessity to take into account a finite size of 
the vortex core, but this certainly demands more numer- 
ical work, while in the present paper the main goal was 
to obtain some analytic results. 
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Appendix A: The calculation of tlie LDOS in the 
nonrelativistic case 



Setting 77 = in Eq. ^ one obtains the solution of 
the Schrodinger equation for B — const without vortex. 
Substituting this solution in the LDOS defintion dTTI) and 
taking into account the widening of the LLs p^ . we 
represent the LDOS as a double sum 



TT ^ — ^ — ^ 

71— Q m— — 00 

1 



En 



(Al) 



where in the second line, we introduced the dimension- 
less variable z — —{E + iT)/Eo with the characteristic 
energy Eq defined below Eq. (fTO|) . To calculate the sum 
in Eq. (jAll) . it is convenient to represent its last factor 
as an exponent 



-S{2n+\m\+7n+l} 



En 



Eoz 



(A2) 



-{l3+S){2n+\m\+m+l)^-l3z 



Here, we also introduced the regularizing exponential fac- 
tor with S > 0, which makes the sum convergent and will 
be set to at the end. Then the LDOS acquires the form 



J 



Im 



dl3e 



00 00 . _2(R-\-S)n 



y' 'e "e 



m— — 00 



Yy[^r'(y)] 



^ r(n+\m\ + 1) 



(A3) 



We operate with the representation (|A3[) in the following we return to the imaginary values z ~{E + ir)/Eo and 

way. First, we consider its analytic continuation for 2: > evaluate the imaginary part. Using Eq. (10.12.20) from 

and perform the calculation. Then to obtain the LDOS, Ref. |32| 

I 



E r(n - 



L^(x)L^(y)z" = (1 - z)-i exp{^z^)ixyz)-^I^ (2^) , \z\ 
■ a + I) I — z \ 1 — z J 



z\ < 1, 



(A4) 



where la is the modified Bessel function, we find the sum 
over n in Eq. (jA3p 
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Nlir,E,B) 



M 



Im 



-ycoth(5+/3) 



m— — oo 



sinh((5 + /3) 



The remaining summation over m in Eq. (IA5|) can 
be done by using the property of the modified Bessel 
function Im{x) = I-rn{x), and that its generating func- 
tion ia^^ 



z™/™(a;)=exp(|[z+l/z]). (A6) 



m— — OO 



We obtain 



Nl{E,B\ 
rini 



(A7) 



Notice that from the last expression one can exphcitly 
observe that it does not depend on y, i.e., in a constant 
magnetic field the LDOS is position independent. Intro- 
ducing a new variable x — 2{5 + /?) we can rewrite the 
last expression as follows: 



Nl(E,B) 



M 



2(7rfi)2 



-Im 



dx- 



2(5 



1 - e-^ 



dx- 



25 



1 - e' 



(A8) 



In the limit (5 — > 0, the second term of Eq. (jASP remains 
real irrespectively the value of z, while the first term gives 



J 



(A5) 



the integral representation of the digamma function: 



38 



^'(z) = -7+ / dt 



1 ~e' 



Re(z) > 0, (A9) 



where 7 is the Euler-Mascheroni constant. Thus we ob- 
tain 



Nl{E,B) 



M 



2{-Kh)' 



-Im 



z + 1 



(AlO) 



so that the final expression for the LDOS after the ana- 
lytic continuation z — —2{E + iT)/{hujc) takes the form 
of Eq. p^ . The oscillatory behavior of the LDOS can 
be explicitly extracted from Eq. (jAlOp [or Eq. (|13p ] using 
the relationship 



^(^-z) = ^{z)A h7rcot(7rz). (All) 

z 



Now we generalize these results for the case when the 
vortex is present. Repeating the steps that led us from 
Eq. ([M]) to Eq. ((X5|) . we obtain 



Nl{v,E,B) 



M 



Im 



dl3e 



-{S+f})-l3z[ 



-acoth(5+^) 
I _ g-2(5+/3) 



rn-\-'q \ 



sinh((5 + (3) 



(A12) 



The sum over m in Eq. (jA12[) is calculated in Ap- was considered above, so that we can concentrate on the 

second term. Substituting it in Eq. (|A12p . we arrive at 
Eq. ^ for AN^{r,E,B). 



pendix IbI Using Eq. (|B7|) we obtain 

-(S+P)(m+r,) r I y 

/ . ^ '-\m+r]\ 



sinh(i5 -h 



X / ^^g-acoshw/sinh(5+,3) 



(A13) 



1 -|_ g-(<5+/3+a;) ■ 



The first term on the right-hand side of the last equa- 
tion corresponds to the LDOS without the vortex, which 



Appendix B: The calculation of the sum over the 
azimuthal quantum number 

The sum over the azimuthal quantum number 
00 

S(r;)= e-''('"+'')/|„+,|(:r). (Bl) 
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can be found using the method described in Ref. |39|. 
Using the integral representation of the modified Bessel 
function^ 



1 

27ri 



(B2) 



where C is a complex path beginning at — ztt + oo and 
ending at in + oo, we obtain 



X cosh uj 



(B3) 



Choosing the contour C to lie in such a way that the 
condition Rew > /3 is satisfied, the series can be made 
convergent, so that 



S(7?) 



2m 



due 



X cosh u 



1 - e-('3+") e("-'5) - 1 



(B4) 

Now changing the variable — > — in the second integral 
we can write 



1 
27ri 



c Jc 



1 a;cosha; 

auje 



-(I3+L0)ri 



(B5) 



where C" is the contour symmetric to the contour C with 
respect to the origin of coordinates. Joining the contours 
C and C leads to an w integral of the form 



duj + 


1 du = 


[ 


Ic 


Jc 





CX3 + 27r 



did 



oo-\-iiT 



did (j) dcj, 
(B6) 

where C" is a rectangle with a length larger than 2(3 and 
width 27rj, centered in the origin, traversed anticlockwise. 
Inside the contour C", the integrand has only one pole 
at Wo = — /3, so that this integral does not depend on 
T] and corresponds to S(0). Therefore, we arrive at the 
final representation for the sum (jBl[) : 



E(r7) 
where 



sm nrj 



due' 



-X cosh LO 



1 + e-('3+") 



E(0), 



E(0) = 6^'=°"'''^. 



(B7) 
(B8) 



Finally, note that one can reproduce the value S(0) from 
Eq. (|B8I) using Eq. (IA6|) . which for z — e~^, reduces to 
the sum (|Bip with 77 = 0. 



Appendix C: The calculation of the density of states 
in the nonrelativistic case 



Substituting Eq. (|15p in the definition ([T5)) and inte- 
grating over the spatial coordinates, we obtain 



ANf^{E,B) = -siuTTTy 



MP 
2{TThy 



-Im 



d/3 







dv 



cosh(w/2) cosh(/3 + S - v/2) 1 + e-" 



(CI) 



where we introduced the new variable v = oj + P + S. This 
double integral can be rewritten using the new variables 

t = e^^^ , X = e"" sls follows: 



MPe 

7^ 



X Im 



7o {l+xY{l + te-^^x) 



(C2) 



where the second integral can be calculated using the 
residue theory 



dxx^ ^ 



{l+x)'^{l + te~'^^x) 



(C3) 



sm TTT] 



(1_ 6-25^)2 



Then the remaining integral is expressed via the hyper- 
geometric function 



dt&'^^^ 
{z + 2t]-1) 

1 + z + 277 



(l-e-2*i)2 



1 - e 



-25 



' 2F,ri,i±i;^;e-- 



1 + z 



1, 



1 + z 



2 ' 2 

3-1-z 



77; e 



-25 



(C4) 



Now we use the series representation of hypergeometric 
functions in Eq. (IC4p : 



e-^*-) / 1 + 2 

2i^l 1, ■ 



z + 2r7 + 1 V 2 

„-25(n-|-7,) 

= E 

n=0 



3 + z 



f], e 



-2S 



z + l + 2ri + 2n 

00 roo 



^ />oo 

/ rfa;e-"(2n+2„+.+i) 



(C5) 



dx- 



-x{z+l+2ri) 



1 - e 



-2x 



where the first one in Eq. (IC4[) is recovered for 77 = 0. We 
observe that the presence of finite S > makes the hyper- 
geometric series well defined, but at the end of the calcu- 
lation the limit (5 — can already be taken. Then taking 
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into account the integral representation of the digamma 
function (|A9p [similarly to Eq. (|A8p ] one can express the 
DOS ()C2|) in the following simple form: 



Ml 



, , z + l \ , ( z + l 
tp 1 — h r; - -0 



(C6) 



which after the analytic continuation z —2{E + 
iT)/{fujJc) takes the final form (|17p . 

Appendix D: Calculation of the function I{y,z,ri) 

As in Ref. i we observe that it is simpler to calculate 
integrals with the derivative dl {y , z , if) / dy representing 
the function /(y, z, rf) in the form: 

f^ dI{Q,z,f^) 
I{y.z,r,) = -j^ —^Q—, (Dl) 

where we used that I{oo,z,ri) — 0. The derivative 
dI{Q, z,vi) / dQ contains two terms; 



dI{Q,z,m) dIi{Q,z,ri) dl2{Q,z,T]) 



dQ 



dQ 



where 



dQ 



-Qcoth,9 



(D2) 



dQ 2 7o sinh^/3 

/•OO 

X / (^^;g-Qcoshi^/sinh^g-(r;-l)w 



dh 1 



-/3(2:+7?-l)^ 



-Qcoth/3 



(D3) 



dQ 2Jo ' sinh^^ 



Using the integral re pres entation of the MacDonald func- 
tion K„{x) (see Ref. 133) 



K^ix) ^ - I e-^™^'^"-'''^dw, (D4) 



we obtain 
dh _ 
dQ JO 

and 

dh ^ 
dQ jQ 



p — Q coth (3 

dp^-l3iz+v)___K^_^{Q/smh(3) 



sinh^ /? 



(D5) 



^ — Qcoth/3 



. X,(g/sinh/3). 
smh p 



(D6) 

Now introducing a new variable t via e — t/{l + t), 
we get 



dQ Jq 



-{z+r,)/2^-2Qt 



(D7) 



X K^_^[2Q^t{l + t)] 



and 

^ = -2e-Q / dtt(^+''-i)/2(l + i)-(-+')-i)/2e-2Qt 



X K^[2Q^t{l + t)]. 



(D8) 



To integrate over t in Eqs. (1071) and (jPSp . we use the 
integral (2.16.10.2) from Ref.[3i 



rP-l 



dx- 



■e-P''K^{c^/x-^ +xz) 

,pz/2 



(a; + z)P 

X T^i/2-p../2 (z+/2) W^i/2-p,./2 iz-/2) , 
^± = ^(p ± y/jP—^), 

Re(p + c) > 0, |arg z| < tt, 2Re(/9) > |Re(i/)|, 

(D9) 



where W\,fi{z) is the Whittaker function. To adapt 
Eq. to the form of Eqs. (|D7| and (|D8]), we have 

to set z = 1, differentiate the result over p, and then take 
the limit p ~> c. This gives 



dx-, —— 

(x + iy 



e-'-^K,[c^x{x + l)] = 



(DIO) 



2' V" ■ 27 'V" 27 " "i/2-P,^/2 ^2 



where the function G\_^{Q) is defined as follows: 

GxAQ) = ^w'AQ) + ^w^AQWi^iQ) ^^^^^ 

+ W;:jQ)W^JQ)-W!,^JQ). 
Accordingly, we obtain that 



^ = r r (ll^ILzl 

dQ \ 2 J \ 2 



^ G(i_2_^)/2,(i_^)/2(Q) 

(D12) 



and 

dh 
dQ 



r ('i^^ r ( G(2_,_,)/2,,/2(Q). 



(D13) 



Now using the differential equation 



W^A')+[-\ + ^ + ^^^^]W,,,iz)^0 (D14) 



and the recursion formula 



z^W^Jz)=[X-^)w,A^) 



1 



Wx-ijz) 



(D15) 
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for the Whittaker function;^ one can transform Gx,^{Q) 
to the form 



M^ + (A-l/2)^ 

LrA^^^yj ^2 

(A -1/2)212 



Q2 

^2- (A -1/2)2 



W^a-i,m(Q) 



VKa,p(Q)VKa-i,p(Q) 



(D16) 



2A- 1 
2Q 



W^l, JQ) - (A - 1/2) 



To obtain the function F\^^j_{y) — — dQG\^^{Q), we 
employ the relationships: 



dQ 
Q 



1 



p-y 

Q 



i^A,p(Q)w^p,M(g) = 



(D17) 



J ^WxAQ)WxMQ) = 



which follow from the differential equation (jD14[) for the 
Whittaker function. Then using the recursion formula 
(jPlSp . we arrive at the following result: 



2i/ 



+ 



2 - (A - 1/2)2 



y 



[wlJy) - Wx-iAvWxAv) - Wx^iAyWx+iAy)] 



2A- 1 
2y 



i^wlAy) - Wx+iAy)9xWxAy) + WxAy)d>^Wx+iAy)]- 



(D18) 



The integral of each term in Eq. (jD2[) is expressed via 
F\,fi{y) with the prefactors given by Eqs. (|D12p and 
(ID13I) . so that we arrive at the final expression (1^^ for 
the function /(y, z, rj) which was defined in Eq. (jDl[) . 

To complete our analysis, we consider the asymptotic 
of the Im/(y, z ^ -{E + iT)/Ea,vi) in the limits y ^ 
and y ^ oo. To do this we use the following representa- 
tions of the Whittaker function: 



r(2/.) 



r(i/2 + /i- A) 



0{y 



(D19) 



and 



W^A,p(y)«e-^V[l + 0(l/y)], y 



oo. 



(D20) 



Substituting Eq. (|D19p in Eq. ([24]) and omitting all real 
terms, which will not contribute to Im/, we obtain 



lral{y -> 0, z, 77) 

TT 



-Im 



4(z-2) 



(D21) 



sinTTjy L V 2 / (z — 1)(2; — 3) 
Now using the property of the digamma function 

tP{z) = i;{z + I) - -, 



(D22) 



we arrive at the result 
lnil{y ^ Q,z,'i]) = 



--; Imtp ^ ^ 

sm T:r] 



(D23) 



One can reproduce the same result directly from Eq. ([23 
Indeed, setting y = in Eq. ([23]) we have 



du- 



(D24) 



18 



Now the integral over lo is elementary and after replacing 
2/3 — T' /?, we obtain 

/(y = 0,z,?7) = - / dp- (D25) 

sm TTi] Jq 1 — e p 

Recognizing in this integral the imaginary part of the 
digamma function (IA9p . we again arrive at Eq. (|D23p . 
The next order corrections to Im/(j/ — >■ 0, z, rj) can be 
obtained by expanding the function dl/dQ at Q = 0, in- 
tegrating the result over Q, and using the y — result 
(jD23[) . For y ^ the expansion contains terms ^ J/^^'' 
and y^ with prefactors that make the resulting approxi- 
mate expression for the LDOS divergent at 77 = 0, 1. 

Substituting the asymptotic (|D20p in Eq. ([M)) . we ob- 
tain that for y 00 



Fx Ay) 



e-yy^^~^ 



+ A(A-1) -^^ 



(D26) 



Appendix E: Solution of the Dirac equation 

A positive energy solution of the time-dependent Dirac 
equation has a form ^'(i,r) = exp{—iEt/h)'^{r), where 
the components of a two-component spinor 



*(r,C) 



V'i(r,C) 
ii/'2(r,C) 



(El) 



satisfy the following equations [compare with Ap- 
pendix A of Ref. 



d_ 

dr 



hvf 



pe-^v I- 



r dip 



he 



^2(r,C) -0, 



d iC d eQA^ 



\or r Oip 



he 



(E2) 



The vector potential A^{r) in Eq. (|E2p is given by 
Eq. ([S]) . From now on, we consider the specific case C = 1 
(omitting the label C in the wave functions) and seek for 
a solution of Eq. (jE2l) in the following form: 

V'i(r) = e'(™-i)'^Vi(0, Mr) = e™'^V'2(r). (E3) 

Then the radial components of the spinor t/^i (r) and ^2 (f) 
satisfy the following system of equations 

hvp ( d m + ri9{r — R) r 
d^ ^ r ^ 2P^ 

d m + rjOir — i?) — 1 r 
r 

(E4) 

Introducing the dimensionless variable y = r^/(2P) and 
denoting the p = R'^/{2P), we rewrite the system (IE4[) 
for y £ [0,p] 

My)-,^^Vy{^y + ^. + -.]My), 



■01 (r) 
Mr) 



E-A 

hvp 

^ E + A \dr 



Since there is no Aharonov-Bohm field for y < p, the 
problem in this domain is identical to that of the Ap- 
pendix D in Ref. Ell. For y e [p, oo[, the system ()E4p 
acquires the form 



hvpV^ ^ ( d 

My) = 7^ — x\jvy{ 



m + T] 1 



■02(2/) 



{E-A)r''\dy ■ 2y 

hv F J— f d m + r; — 1 
~{E + A)l^\d:y Yy 



My). (E6a) 



■01 (l/)- 



(E6b) 



The matching conditions (PSI) and take the form 



V'i(p-o) = v^i(p + o), 

^[{p + Q)-^[{p-Q) = ^^^{p) 

2p 



(E7) 



and 



^2(p-0) =7A2(p + 0), 

Mp + ^)-Mp-^)- 



2p 



Mp). 



(E8) 



where the derivative is taken over y. One can obtain from 
the system (|E5|) that for y £ [0, p] the spinor components 
satisfy the following second-order differential equations: 



rf2 Id 1 (m - 1)2 _^ A - 771 

dy2 y dy 4 



2y 



d"^ Id 1 771^ A-771-fl 

dy"^ y dy 4 47/^ 27/ 



My) = Q, 

(E9a) 

^2(2/) = 0, 

(E9b) 



Mr) 



where we introduced A — (E^ — A'^)P /(hvp)"^ ■ The sec- 
ond order differential equations for the domain y £ [p,oo[ 
corresponding to the system (IE6P can be obtained from 
Eq. (|E9p by replacing m ^ m + rj: 



My)^o, 

(ElOa) 

My) = 0- 

(ElOb) 



Id 1 (777 + 7^—1)2 \ — m — rj 
dy'^ ij dy 4 



Mr). 



rf2 1 d 

dy^ y dy 4 



4y2 2y 

1 (777 + 7;)2 A — 777 — 77 -|- 1 



47/2 



22/ 



(i; - A)l 

hvFV2 ^ f d 



777 1' 
777 —1 



y 



The equations can be reduced to the equations for the 
degenerate hypergeometric function (see Eq. (6.3.1) of 
Ref. Hi) and the solutions of Eqs. (|E9al) and (lElOap are 
given, respectively, by 



i^i{y)- 



(E5) 
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m— ll+m + l — A 



1 + |m - l|;y , r <R 



My) =A,r.y^'^+^~^\'\ryl^^ ( ^+-^, 1 + |m + 7? - 1|; y 



- A 



1 - |m + 77 - l|;y , r > R, 



(Ella) 



(Ellb) 



where a± = rn + 77 + 1 ± |m + ?; — 1|, A,n, and C,„ 
are constants, and $(a,c;z) is the confluent hypergeo- 
metric function. The solution (jEllap contains only one 
term due to the condition of square integrability and the 
absence of the Aharonov-Bohm field for r < R. Writing 
the solution (|Ellb|) we used that for noninteger c the so- 
lution of Eq. (jE10a[) can be expressed via ^{a,c;z) and 
2:i-c$(a _ c + 1, 2 — c; z). The coefficients Am, Bm, and 
Cm to be found from the matching conditions (jETj. The 
consideration of the limit i? — (p — 0) greatly simpli- 
fies the calculation, because one can expand the solutions 
to the linear in p terms. Then one finds that 



My) - Amy^-'+^'-'^/'e-y/' 

,'|m + 7?— l|+m + r? + l — A 
X $ ' 



l + |r7i-)-77— 1|;?; 

(E12) 



Since ^{a,c;z) behaves as at large y unless a = —n 
with n = 0,1,2,..., in order to have the square inte- 
grable solutions, the value A should be equal to the eigen- 
value Xm.n defined by Eq. ([33)1 . In this case $ function 
is reduced to the generalized Laguerre polynomials [see 
Eq. (6.9.2.36) of Ref.lHl 

Kiy) = ^^^^^^-n,a + l,y). (E13) 
L [a + l)n] 

Introducing the functions 

1 /2 

one can rewrite the solution (jE12l) in a more compact 
form 



The definition (jE14|) generalizes the functions considered 
in Ref. 42 for the case of the noninteger i' > —1. These 
functions satisfy the following orthogonality condition 



(E16) 



Then having ipi{y) one can find ip2{y) from Eq. (|E6bp 
using the recursion formulas^ 



(a; + z.)j;(x) = 

2x{dldx),r^{x) = 
[x(n + i/)]5j^_^(a:) 



[x{7i + v + l)]^Jl\^{x), 



[x{n + v + l)]^-J1^M^)- 



(E17) 



Then demanding that the spinors obey the normalization 
condition 



27r 



My) = AmJlm+r^_i\{y)- 



(E15) 



dip I rdr*;,„,(r, C)*nm(r, C) = ^n,n"5m.m' 

Jo 

(E18) 

we obtain the solutions (PO]) . (l5Tt . and ([5^ for n > 0. 
The zero-mode solutions have to be considered sepa- 
rately. Analyzing the initial system (jE6|) . one finds that 
the only allowed solution is the negative energy i? = — A, 
TO < with ipi(r) = 0. The corresponding spinor is given 
by Eq. (p4)) . One can verify that for 77 = these solu- 
tions transform up to the phase factors to the solutions 
obtained in Ref. 41. To show this, one should relabel the 
quantum numbers n + m ^ n for m > 1 and n -f 1 — > n 
for TO < and use the property J"(y) = (— 1)" J"J'^(y), 
which is valid only when J^jy) is defined for the inte- 
ger values of v as done in Ref. [i^. After this relabeling 
is made, the spectrum (j33p acquires a conventional form 
dependent only on the LL index which for A = reduces 
to Eq. ©. 
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